; analytic xHI solution for tau0 * xHI >> 1
;-------------------------------------------------------------

function x_tau0_x_big, nH, T, GHI, tau0, y, dbg=dbg

  if n_elements(dbg) eq 0 then dbg = 0

  nH = double(nH)
  T  = double(T)
  GHI = double(GHI)
  tau0 = double(tau0)
  y = double(y)

  CI = hici_hui( T )
  RC = hiirca_hui( T )

  CI = double(CI)
  RC = double(RC)

  RR = -( GHI * tau0 * exp(-tau0) + (CI + RC) * nH )
  QQ = GHI * exp(-tau0) * (tau0-1.0d0) + CI * nH - (CI+RC) * nH * y 
  PP = GHI * exp(-tau0) + CI * nH * y
  d = QQ*QQ - 4 * PP * RR

  q = - 0.5d0 * (QQ + abs(QQ)/QQ * sqrt(d))
  
  xHIIp = q / RR      ; this one
  xHIIm = PP / q                

  xHIp = 1.0d0 - xHIIp

  return, xHIp

end
